library(parallel)
library(doMC)
library(DHARMa)
library(broom.mixed)
library(MuMIn)
library(mgcv)
library(bbmle)
library(lme4)
library(merTools)
library(magrittr)
library(gridExtra)
library(ggplot2)
library(tidyr)
library(purrr)
library(plyr)
library(dplyr)

1 - Variables

table 1: Variáveis utilizadas na descrição estatística

Variables
code description class range
n_nRef number of SADs not refuted by the KS bootstrap test with critical p-value 0.05 continuous (0 ; 100)
diffS_mean diffS = (S_MN - S_obs) / S_obs, diffS_mean = mean(diffS); S_MN = species richness estimated by MN continuous (-0.977 ; 4.78)
U_med average of 10 replicates of the speciation rate estimated by MNEE continuous (8.860e-5 ; 4.221e-2)
p proportion of tree cover continuous (0.013 ; 0.961)
MN Neutral Model (EE = spatial explicit; EI = spatial implict) category 2 levels
d mean dispersal distance (meters) continuous (0.456 ; 19.183)
k proportion of propagules until the nearest neighborhood seq(0.99 : 0.05) category 20 levels
d_Lplot mean dispersal distance / side of the sample area continuous (0.02 ; 0.192)
S_obs observed species richness integer (26 ; 458)
Ntotal number of individuals in the sample area integer (649 ; 12105)
SiteCode sample area code category 103 levels

Figure 1 Possible Predictor Variables

2 - Number of unrefuted SADs

Figura 2.1 número de SADs não refutadas ~ p * k * MN. A linha azul é uma estimativa baseada em ‘loess’.

Figura 2.2 número de SADs não refutadas ~ (d * MN|Site ~ p). Site está ordenado pelo valor de p.

2.1 GLMM binomial

Modelo cheio 1:

linear term

  1. ~ (p + p^2) * (d + d^2) * MN; if : var_dispersao = contiguous
  2. ~ ( p + p^2) * k * MN; if : var_dispersao = category

random term

  1. 1|SiteCode
  2. MN|SiteCode
  3. (var_dispersao + var_dispersao^2)*MN|SiteCode, if : var = dispersao : contiguous
# dados
df_md <- df_resultados %>% 
  select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>% 
  gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z) 
# funcao ajuste dos modelos
f_md <- function(dados){
  var_dispersao <- unique(dados$var_dispersao)
  if(var_dispersao == "k"){
    l_md <- vector("list",2)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"))
    dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
  }else{
    dados$value_dispersao <- as.numeric(dados$value_dispersao)
    l_md <- vector("list",3)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"),
                     paste0(var_dispersao," (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         ( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    
  }
  return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio) <- NULL
l_nRef__modeloCheio <- do.call(c,l_nRef__modeloCheio)
#
f_warning <- function(glmm_object){
  # glmm_object <- l_nRef__modeloCheio[[4]]
  v_message <- glmm_object@optinfo$conv$lme4$messages %>% 
    as.character()
  if(length(v_message)==0){v_message <- "OK"}
  if(length(v_message)>1){v_message <- v_message[1] }
  return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio,f_warning,.id="glmer") %>% 
  rename(warning_message = V1)

Table 2.1 Modelos Cheios estimados e avisos de convergência

glmer warning_message
d_Lplot.z 1|Site Model failed to converge with max|grad| = 0.00343548 (tol = 0.002, component 1)
d_Lplot.z MN|Site Model failed to converge with max|grad| = 0.00836802 (tol = 0.002, component 1)
d_Lplot.z (d_Lplot.z+d_Lplot.z^2)*MN|Site OK
d.z 1|Site Model failed to converge with max|grad| = 40.1827 (tol = 0.002, component 1)
d.z MN|Site Model failed to converge with max|grad| = 0.00746139 (tol = 0.002, component 1)
d.z (d.z+d.z^2)*MN|Site OK
k 1|Site OK
k MN|Site OK

allFit

i <- 1
names(l_nRef__modeloCheio[v_glmerUpdate])[i]
#
## update1:
md_allFit <- allFit(l_nRef__modeloCheio[v_glmerUpdate][[i]],maxfun = 1e5, parallel = 'multicore', ncpus = 3)

1) d_Lplot 1|Site

7 optimizer(s) failed

2) d_Lplot MN|Site

7 optimizer(s) failed

3) d 1|Site

7 optimizer(s) failed

4) d MN|Site

7 optimizer(s) failed

Modelo cheio 2:

# dados
df_md <- df_resultados %>% 
  select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>% 
  gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z) 
# funcao ajuste dos modelos
f_md <- function(dados){
  var_dispersao <- unique(dados$var_dispersao)
  if(var_dispersao == "k"){
    l_md <- vector("list",2)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"))
    dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN +
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
  }else{
    dados$value_dispersao <- as.numeric(dados$value_dispersao)
    l_md <- vector("list",4)
    names(l_md) <- c(paste0(var_dispersao," 1|Site"),
                     paste0(var_dispersao," MN|Site"),
                     paste0(var_dispersao," ",var_dispersao,"*MN|Site"),
                     paste0(var_dispersao,"^2 (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
    l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN + 
                         (1|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN + 
                         (MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * value_dispersao * MN + 
                         (value_dispersao * MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    l_md[[4]] <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                         (p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN + 
                         ( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
                       family = "binomial",data=dados,
                       control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
    
  }
  return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio2 <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio2) <- NULL
l_nRef__modeloCheio2 <- do.call(c,l_nRef__modeloCheio2)
#
f_warning <- function(glmm_object){
  v_message <- glmm_object@optinfo$conv$lme4$messages %>% 
    as.character()
  if(length(v_message)==0){v_message <- "OK"}
  if(length(v_message)>1){v_message <- v_message[1] }
  return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio2,f_warning,.id="glmer") %>% 
  rename(warning_message = V1)

Table 2.2 Modelos Cheios estimados e avisos de convergência

glmer warning_message
d_Lplot.z 1|Site Model failed to converge with max|grad| = 0.00284637 (tol = 0.002, component 1)
d_Lplot.z MN|Site OK
d_Lplot.z d_Lplot.z*MN|Site OK
d_Lplot.z^2 (d_Lplot.z+d_Lplot.z^2)*MN|Site OK
d.z 1|Site Model failed to converge with max|grad| = 0.00369235 (tol = 0.002, component 1)
d.z MN|Site OK
d.z d.z*MN|Site OK
d.z^2 (d.z+d.z^2)*MN|Site OK
k 1|Site OK
k MN|Site OK

allFit

d_Lplot.z 1|Site
7 optimizer(s) failed

d.z 1|Site
7 optimizer(s) failed

Comparação de Modelos Cheios:

Tabela 2.3 Comparação baseada em AICc dos modelos cheios

GLMM dAICc df weight
d^2 (d+d^2)*MN|Site 0.0 39 1
d_Lplot^2 (d_Lplot+d_Lplot^2)*MN|Site 274.6 39 <0.001
d d*MN|Site 38427.3 22 <0.001
d_Lplot d_Lplot*MN|Site 38586.5 22 <0.001
k MN|Site 74141.2 123 <0.001
d MN|Site 95710.0 15 <0.001
d_Lplot MN|Site 96279.6 15 <0.001
k 1|Site 106158.0 121 <0.001

Tabela 2.4 Coeficiente de Determinação Condicional e Marginal

##                   R2m       R2c
## theoretical 0.3232656 0.9128366
## delta       0.3111874 0.8787302

Diagnostico do modelo cheio mais plausível

Figura 2.3 Resíduos Quantílicos do modelo cheio plausível

Figura 2.4 Quantile-quantile plot random effects.

Figura 2.5 Predito e observado pelo modelo. Em vermelho linha y=x; em azul um modelo linear por sítio de amostragem.

Figura 2.6 Predito e observado para modelo cheio

Seleção da combinação de preditoras:

Modelo Global

md_nRef <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                   (p.z + I(p.z^2)) * (d.z + I(d.z^2)) * MN + 
                   ( (d.z + I(d.z^2) ) * MN|SiteCode),
                 family = "binomial",data=df_resultados,
                 control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)),
                 na.action = "na.fail")

Tabela 2.5 Tabela de Seleção de Modelos, deltaAICc<7

## Global model call: glmer(formula = cbind(n_nRef, 100 - n_nRef) ~ (p.z + I(p.z^2)) * 
##     (d.z + I(d.z^2)) * MN + ((d.z + I(d.z^2)) * MN | SiteCode), 
##     data = df_resultados, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
##         optCtrl = list(maxfun = 2e+09)), na.action = "na.fail")
## ---
## Model selection table 
##        (Int)     d.z   d.z^2 MN    p.z   p.z^2 d.z:MN d.z:p.z d.z:I(p.z^2)
## 35712  2.296 -0.3911 0.03539  + 0.4903 -0.3671      +  0.5471             
## 35840  2.382 -0.2154 0.03597  + 0.4569 -0.4532      +  0.4794     -0.17770
## 35696  1.929 -0.3912 0.03512  + 0.6350              +  0.5464             
## 43904  2.296 -0.3940 0.03634  + 0.4741 -0.3668      +  0.6528             
## 44032  2.383 -0.2168 0.03648  + 0.4409 -0.4545      +  0.5853     -0.17890
## 43888  1.929 -0.3939 0.03581  + 0.6170              +  0.6526             
## 36736  2.265 -0.3906 0.07141  + 0.5017 -0.3356      +  0.5468             
## 39808  2.435 -0.3912 0.03544  + 0.4367 -0.5064      +  0.5470             
## 44928  2.264 -0.3932 0.07143  + 0.4864 -0.3345      +  0.6518             
## 36864  2.363 -0.2319 0.04923  + 0.4622 -0.4336      +  0.4857     -0.16090
## 39936  2.411 -0.2225 0.03575  + 0.4459 -0.4824      +  0.4821     -0.17040
## 48000  2.436 -0.3939 0.03617  + 0.4201 -0.5077      +  0.6526             
## 45056  2.364 -0.2318 0.04964  + 0.4481 -0.4352      +  0.5912     -0.16350
## 48128  2.414 -0.2237 0.03674  + 0.4291 -0.4855      +  0.5885     -0.17220
## 40832  2.380 -0.3905 0.06743  + 0.4577 -0.4514      +  0.5465             
## 106368 2.501 -0.3906 0.08391  + 0.4095 -0.5720      +  0.5464             
## 49024  2.381 -0.3932 0.06788  + 0.4410 -0.4525      +  0.6519             
## 56320  2.389 -0.2016 0.03603  + 0.4525 -0.4603      +  0.4841     -0.19170
## 40960  2.392 -0.2376 0.04888  + 0.4503 -0.4610      +  0.4868     -0.15380
## 114560 2.502 -0.3933 0.08434  + 0.3934 -0.5732      +  0.6523             
## 64512  2.423 -0.2323 0.03672  + 0.4255 -0.4949      +  0.5918     -0.16360
## 106496 2.478 -0.3213 0.07180  + 0.4188 -0.5488      +  0.5200     -0.07021
## 57344  2.384 -0.2154 0.04265  + 0.4555 -0.4550      +  0.4866     -0.17740
## 114688 2.479 -0.3211 0.07214  + 0.4020 -0.5501      +  0.6254     -0.07340
## 65536  2.408 -0.2846 0.06205  + 0.4309 -0.4792      +  0.6110     -0.11000
##        I(d.z^2):MN I(d.z^2):p.z d.z^2):I(p.z^2 MN:p.z I(p.z^2):MN
## 35712            +     -0.04297                     +            
## 35840            +     -0.04304                     +            
## 35696            +     -0.04152                     +            
## 43904            +     -0.07235                     +            
## 44032            +     -0.07353                     +            
## 43888            +     -0.07158                     +            
## 36736            +     -0.05651      -0.036240      +            
## 39808            +     -0.04324                     +           +
## 44928            +     -0.08643      -0.035770      +            
## 36864            +     -0.04843      -0.013710      +            
## 39936            +     -0.04349                     +           +
## 48000            +     -0.07313                     +           +
## 45056            +     -0.07795      -0.013140      +            
## 48128            +     -0.07324                     +           +
## 40832            +     -0.05550      -0.032450      +           +
## 106368           +     -0.06169      -0.048910      +           +
## 49024            +     -0.08515      -0.032170      +           +
## 56320            +     -0.04582                     +           +
## 40960            +     -0.04801      -0.013600      +           +
## 114560           +     -0.09139      -0.048620      +           +
## 64512            +     -0.07326                     +           +
## 106496           +     -0.05709      -0.036540      +           +
## 57344            +     -0.04806      -0.007001      +           +
## 114688           +     -0.08638      -0.035920      +           +
## 65536            +     -0.08312      -0.025920      +           +
##        d.z:MN:p.z d.z:I(p.z^2):MN I(d.z^2):MN:p.z I(d.z^2):I(p.z^2):MN df
## 35712                                           +                      33
## 35840                                           +                      34
## 35696                                           +                      32
## 43904           +                               +                      34
## 44032           +                               +                      35
## 43888           +                               +                      33
## 36736                                           +                      34
## 39808                                           +                      34
## 44928           +                               +                      35
## 36864                                           +                      35
## 39936                                           +                      35
## 48000           +                               +                      35
## 45056           +                               +                      36
## 48128           +                               +                      36
## 40832                                           +                      35
## 106368                                          +                    + 36
## 49024           +                               +                      36
## 56320                           +               +                      36
## 40960                                           +                      36
## 114560          +                               +                    + 37
## 64512           +               +               +                      37
## 106496                                          +                    + 37
## 57344                           +               +                      37
## 114688          +                               +                    + 38
## 65536           +               +               +                      38
##           logLik    AICc delta weight
## 35712  -14179.02 28424.6  0.00  0.113
## 35840  -14178.06 28424.7  0.11  0.107
## 35696  -14180.31 28425.1  0.54  0.086
## 43904  -14178.29 28425.2  0.58  0.085
## 44032  -14177.32 28425.3  0.66  0.081
## 43888  -14179.58 28425.7  1.13  0.064
## 36736  -14178.68 28425.9  1.34  0.058
## 39808  -14178.82 28426.2  1.64  0.050
## 44928  -14177.96 28426.5  1.94  0.043
## 36864  -14178.03 28426.7  2.08  0.040
## 39936  -14178.05 28426.7  2.13  0.039
## 48000  -14178.09 28426.8  2.21  0.037
## 45056  -14177.28 28427.2  2.62  0.031
## 48128  -14177.31 28427.3  2.67  0.030
## 40832  -14178.55 28427.7  3.13  0.024
## 106368 -14177.79 28428.2  3.64  0.018
## 49024  -14177.83 28428.3  3.71  0.018
## 56320  -14178.00 28428.7  4.06  0.015
## 40960  -14178.04 28428.7  4.14  0.014
## 114560 -14177.06 28428.8  4.22  0.014
## 64512  -14177.30 28429.3  4.69  0.011
## 106496 -14177.72 28430.1  5.53  0.007
## 57344  -14177.99 28430.7  6.08  0.005
## 114688 -14177.01 28430.8  6.16  0.005
## 65536  -14177.20 28431.1  6.55  0.004
## Models ranked by AICc(x) 
## Random terms (all models): 
## '(d.z + I(d.z^2)) * MN | SiteCode'

Model Averaging: